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In  the  context  of  a  model  problem  we  describe  post-processing  techniques 
for  the  calculation  of  generalized  stress  intensity  factors.  We  discuss  two 
broad  classes  of  methods,  one  involving  an  "influence”  function,  and  the 
other  related  to  the  well-known  energy  release  principle  of  fracture  mechanics. 
An  error  analysis  is  sketched  and  two  numerical  examples  are  given  to  illustrate 
the  effectivity  of  the  techniques. 
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Abstract 

T^L  Cudlr*' 

_ >In  the  context  of  a  model  problem  we"  describe  post-processing 

techniques  for  the  calculation  of  generalized  stress  intensity 
factors.  ^Wt^discuss  two  broad  classes  of  methods,  one  involving 
an  influence^  function,  and  the  other  related  to  the  well-known 
energy  release  principle  of  fracture  mechanics.  An  error  analysis 
is  sketched  and  two  numerical  examples  are  give  to  illustrate  the 
effectivity  of  the  techniques. 
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51.  Introduction 

51.1.  General  Introduction 

This  is  the  second  in  a  series  of  three  papers  in  which  we 
discuss  "post-processing"  as  it  may  be  applied  in  the  finite  element 
method.  In  the  first  paper  [  1  ]  we  briefly  described  some  general 
aspects  of  post-processing,  and  illustrated  a  few  of  the  ideas  in 
the  particular  cases  of  two  model  problems  from  structural  mechan¬ 
ics.  For  these  model  problems  we  were  concerned  with  obtaining 
values  for  the  displacements,  stresses,  bending  moments  etc., 
either  at  a  point  or  as  an  average  over  some  subsection  of  the 
structure.  However,  we  were  careful  to  note  that  for  two  dimen¬ 
sional  problems  the  methods  discussed  could  only  be  expected  to 
perform  well  if  the  point  or  subsection  under  consideration  was 
"reasonably  distant"  from  any  corner  point  of  the  region  occupied 
by  the  structure.  In  this  paper  we  take  up  the  problem  of  how  to 
proceed  when  this  "reasonably  distant"  criterion  is  no  longer 
satisfied. 

As  is  well  known,  corner  points  usually  give  rise  to  some 
form  of  singular  behavior  in  the  derivatives  of  the  displacements. 

In  fact,  the  stress  intensity  factors  of  elasticity  theory  are 
the  coefficients  of  the  terms  that  exhibit  this  singular  behavior. 
They  serve  as  a  means  of  characterizing  the  state  of  stress  in 
the  neighborhood  of  the  corner  point.  We  shall  see  that  calcu¬ 
lation  of  these  stress  intensity  factors,  as  well  as  certain 
coefficients  of  non-singular  terms,  is  possible  within  the  post¬ 
processing  scheme.  If  desired,  values  of  the  displacements  or 
stresses  at  points  near  the  corner  may  then  be  found  by  using 
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the  values  of  the  stress  intensity  factors, along  with  the  non¬ 
singular  coefficients, in  known  asymptotic  expansions  valid  in 
the  vicinity  of  the  corner  point. 

In  this  paper  our  post-processing  calculations  will  be  based 
on  expressions  that  are  of  a  slightly  more  general  form  than  those 
dealt  with  in  the  first  paper  of  this  series.  In  place  of  the 
model  expression  for  the  generic  quantity  4>  which  appears  there 
(see  (1.1)  of  [1  ]),  we  will  now  consider  expressions  for  4» 
which  take  the  form 

(1.1)  $  =  *(u>)  =  f  X  •  VwdA  +  f  £u>dA  +  f  pa)ds  +  R, 

Ja  h  ha 

where  X  =  (X^,!^)  is  a  vector  function  on  (2,  C  is  a  scalar 

function  on  ft,  p  is  a  scalar  function  defined  on  3(2,  or  some 

specified  portion  of  3(2,  and  R  is  an  integral  which  is  easily 
computable  from  the  load  data  of  the  problem.  If  w  is  a  finite 
element  approximation  to  uj,  then  (1.1)  suggests  that  we  use 

(1.2)  $  =  $(w)  =  [  X  •  VudA  +  f  £udA  +  [  puds  ♦  R 

Ja  1 (2  '  3(2 

as  an  approximation  to  $  (c.f.  (1.2)  of  Cl]). 

§1.2.  formulation j^fjthe  n^del^groblem 

As  usual,  we  shall  describe  the  post-processing  technique  in 
the  context  of  a  particular  model  problem.  Suppose  (2  is  a 
bounded,  two  dimensional  region  whose  boundary  is  made  up  of  smooth 
arcs  r.,...,r  .  Consider  a  boundary  value  problem  governed 

by  the  equation 
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(1.3a)  V2u>  =  0  in  0, 

and  for  which,  independently  on  each  I*.,  either  the  Dirichlet 
boundary  condition 

(1.3b)  <i>  =  f.  on  r. 

J  3 

or  the  Neumann  boundary  condition 

A  A 

(1.3c)  (Vw)  •  n  =  g.  on  I*.  (n  =  outward  pointing  unit  normal) 

applies.  For  definiteness,  assume  that  and  ^  are  adjacent 

edges  of  0  which  meet  at  the  origin  of  the  coordinate  system 
(x^,x2).  Establish  polar  coordinates  (r,0)  at  this  origin,  and 
suppose  that  and  are  in  fact  straight  line  segments 

which  correspond  to  6=0  and  6  =  an  respectively.  Assume 
that  ,in  the  vicinity  of  0  the  region  J)  corresponds  to  the  cone 
0  <  0  <  an  (see  Figure  1).  For  and  let  the  boundary 

conditions  (1.3b)  and  (1.3c)  take  the  specific  form 

(1.3d)  (o  *  f .  *  0  on  r^»  and  ( Vw)  •  n  =  g^  =  0  on  r2* 

(If  a  =  2,  then  what  we  have  described  models  a  slit  domain. 

The  upper  face  (6=0)  of  the  slit  is  fixed,  whereas  the  lower 
face  (0  =  2n)  is  traction  free.) 


Figure  1. 


In  a  neighborhood  of  the  vertix  0  the  solution  o>  and  its 
derivatives  have  asymptotic  expansions  of  the  form 


(1.4a) 


u> 


N 

„  h 


( 2n-l) 


k  r 
n 


7Z 


sin( ( 2n-l)  ^ ) 


(2N+D^ 
0(r  2ot), 


N 


5n  ,  (2n-l-2a)i/sin((2n-1-2“,^)' 

(1.4b)  Vw  =  l  k  r  2a‘ 

n=l  2a  n 


/sin((2n-l-2a)^)\  ,  ( 

r  °(r 

''Jos((2n-l-2a)^-)/  ' 

\  2a  / 


2N+l-2a) 


for  some  numbers  k  ,n  =  l,...,N,  with  N  arbitrary.  Provided 
the  fj  and  gj’s  are  sufficiently  smooth,  then  in  the  vicinity 
of  any  other  vertex,  P  say,  of  fl 
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uj  =  0(1) 

(1.5) 

va >  =  o(  |x-x(p)  r3/4) 

where  |x-x( P)  |  =  [  (x^-x^P))  2  +  (x1~x2 (P)  )  2]1/ 2 .  See  [  2  ]  for 
further  details.  Notice  that  if  n  <  y  +  a,  the  n**1  term  in 
the  expansion  (1.4b)  for  is  unbounded  as  r  0.  We  can 

think  of  the  coefficients  k^  of  these  terms  as  analogs  of  the 
stress  intensity  factors  of  elasticity.  In  this  paper  we  general¬ 
ize  this  terminology  somewhat,  and  refer  to  all  the  coefficients 
kn,  whether  or  not  the  corresponding  terms  in  (1.4)  are  singular, 
as  stress  intensity  factors. 

Our  approach  to  obtaining  approximate  values  for  the  "dis¬ 
placement"  uj  and  the  "stress"  near  0  will  be  to  first 

o  Xj 

calculate  the  coefficients  k^  b.y  a  post-processing  of  the 
finite  element  solution,  and  then  use  the  asymptotic  expansions  (1.4) 
to  find  the  displacement  and  stresses.  Notice  that  this  is  a 
reversal  of  a  "direct"  method  which  is  often  used  to  calculate 
the  stress  intensity  factors:  firstly,  the  finite  element  solution 
is  "fitted"  in  some  way  to  (1.4),  and  the  resulting  coefficients 
corresponding  to  the  kn  are  taken  as  the  appropriate  stress 
intensity  factors. 

Although  our  considerations  in  this  paper  will  be  restricted 
to  the  above  model  problem  the  ideas  of  our  analysis  extend  in  a 
natural  way  to  many  other  two  dimensional  problems.  In  particular, 
everything  we  do  could  also  be  carried  out  in  the  context  of  two 
dimensional  elasticity. 
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i  1 . 3 .  OutJLine, jof^the. jfiager 

In  §2  and  §3  we  derive  two  different  expressions  for  the 
stress  intensity  factors  k^.  These  expressions  are  both  of 
the  general  form  described  by  (1.2).  Section  4  addresses  the 
accuracy  of  the  post-processed  approximations  to  k^  based  on 
these  expressions.  Our  discussion  at  this  point  will  be  very 
similar  to  that  in  the  corresponding  section  of  our  earlier 
paper.  Finally,  in  §5  we  give  two  numerical  examples  to 
illustrate  the  practical  effectivity  of  post-processing  in  the 
current  context . 
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§2  .  Thje  Jlerie^a f  Iai ejn ce  c t  ^or^J^e t hod 

§2.1. 

The  first  method  for  calculating  stress  intensity  factors 
that  we  shall  describe  is  based  on  an  influence  function  approach . 
At  least  in  theory,  the  concept  of  an  influence  function  for  a 
stress  intensity  factor  has  been  known  for  some  time  (see,  for 
instance,  [  3]).  However,  these  kinds  of  ideas  do  not  appear  to 
be  widely  used  ’n  finite  element  practice  (though,  see  [4  ],[5]  and 
[10]).  Although  what  we  shall  say  will  only  be  in  the  context 
of  the  model  problem  (1.3),  the  ideas  may  be  readily  extended  to 
other  situations. 


§2.2.  ^njexgress^ion  ^for  km 

Let  <p  be  a  function  defined  on  £2.  For  the  moment,  let  us 

only  assume  that  <p  is  smooth  everywhere  in  £2  except,  perhaps, 

near  the  corner  of  interest  at  0.  Let  ( £  =  0,...,s-l)  be 

some  enumeration  of  the  vertices  of  £2,  with  denoting  the 

vertex  0.  Fox-  e  >  0,  sufficiently  small,  remove  from  £2  discs 

of  radius  e  about  each  P^ .  Denote  by  £2e  the  new  region  so 

formed.  Let  and  Y#e^  be  as  indicated  in  Figure  2.  Multiply 

J  * 

(1.3  a)  by  <p  and  integrate  by  parts  over  £2^  to  obtain 


(2.1)  0=[  (V^uj)<pdA  =  (T  [  +  7  [  )  (Vwn<p-V(p*  nuOds 

JfL  *  J  Ce)  j  J  (e> 

e  Y*  j 

+  f  (V2<p)wdA, 

J  n 


A 

where  n  is  the  outward  pointing  unit  normal. 
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Figure  2. 


Now,  let  us  place  some  further  restrictions  on  cp.  Suppose 

that 

(I)  <p  a  o  on  each  of  the  edges  Tj  to  which  the  Dirichlet 
boundary  condition  (1.3b)  applies. 

(II)  In  the  vicinity  of  the  vertex  0, 


-(2m-l)^ 


it  (  2m- 1  ] 


sin((  2m-l)  ttt)  +  <p. 


for  some  m  =  1,2...  and  some  smooth  function  <p^  (bounded 

second  derivatives  will  suffice) . 

Using  the  above  properties  of  cp,  along  with  the  expansions 

(1.4)  and  (1.5),  we  may  directly  evaluate  the  line  integrals  over 
( e) 

the  y  ’  s  appearing  in  (2.1).  This  gives 
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is  bounded  near  0, 

leads 

to  the  following 

expression  for  k^: 

f  0 

N  / 

r 

D  r 

(2,2)  k  =  (V  (p)u>c 

m  _ 

l 

j  J 

(g  *  <p-V<p*  na))ds  - 
ri 

I  7(p*nf.ds, 

3  f-  3 

where 


N  D 

£  and  £  denote  summation  over  all  the  edges 
3  3 


for 


which  the  Neumann  boundary  condition  (1.3c)  and  the  Dirichlet 
boundary  condition  (1.3b)  apply  respectively.  Notice  that  if 
<P  had  the  additional  properties 


(i)  V2<P 


0  in  ft ,  and 


(ii)  V<p  •  n  =  0  on  the  edges  I*, 
condition  (1.3c)  applies, 


for  which  the  Neumann 


then  (2.2)  would  become  simply 


(2.3) 


N  , 

l  g^ds  - 

i  Jr.  J 

J  3 


V<p  •  nf .  ds  . 
3 


The  function  <p  could  then  be  thought  of  as  an  influence  function, 

relating  the  applied  boundary  displacement  and  traction  loading 

to  the  resulting  stress  intensity  factor  k  at  0  (see  [33). 

m 

In  general  such  a  <p  is  not  immediately  available. 

The  expression  (2.2)  is  of  the  form  (1.1),  and  so  in  line 

with  (1.2)  we  consider  the  approximation  k  to  k  calculated 

mm 

using  the  finite  element  solution  oj, 


(2.4) 


'N/ 


f  2  ^  f  u  f 

(V  cp)  (i)dA+E  ( g^P  ~V<p  •  nus)ds  -  £  |  V<p  *nf  .ds . 

)«  j  3  j  Jfj  ^ 


Much  as  in  [  1  ],  we  shall  see  in  §4  that  the  accuracy  of  km 

2  A 

is  influenced  by  the  smoothness  of  £  =  V  <p  and  p=  -V<p  •  n. 

As  before,  the  issue  of  how  to  select  a  <j>  for  use  in  (2.4) 
is  essentially  a  matter  of  deciding  how  to  satisfy  the  boundary 
requirement  (I)  above  in  a  smooth  way.  Among  others,  the  cut 
off  function  and  blending  function  techniques  described  in  [  1  ] 
may  be  used.  We  shall  leave  to  §5  any  further  discussion  of 
the  actual  numerical  implementation  of  (2.4). 


§ 3 .  The  Generalized  Energy  Release  Method 
5  3.1.  Pri el^minaries 

The  second  approach  that  we  shall  discuss  is  related  to  the 
well-known  energy  release  method  of  fracture  mechanics.  Over  the 
years  this  method  has  been  implemented  in  numerous  forms.  To  name 
but  a  few:  the  J-integral  [6,9,10],  the  stiffness  derivative  [7] 
and  the  crack  closure  [8,9,10]  methods.  In  common  with  these 
methods,  the  approach  that  we  shall  outline  in  this  section  only 
applies  to  cases  where  the  included  angle  at  the  vertex  0  is 
it  or  2 tt  (i.e.,a=l,2).  This  restriction  is  not  too  serious 
since  in  practice  these  cases  are  by  far  the  most  important. 
Again,  though  we  only  discuss  the  model  problem  (1.3),  the 
technique  that  we  shall  outline  may  be  extended  to  treat  other 
problems . 


§3.2.  Another,  jej^res^sicsn^fj or  k 

Let  cp  be  a  smooth  function  on  ft  (bounded  first  derivatives 
on  n  will  suffice),  and  define 


(3.1)  v  = 


_ -4a _ 

*  (  2m-l)(l+ 2a-2m)  r 


( l+2a-2m)-^~ 
2a 


sin((l+2a-2m)^-  )  (m  =  1,2, 


For  e  >  0,  sufficiently  small,  and  with  fi  ,  and  Yp£^ 

e  ] 


as  in  52.2,  consider  the  integral 
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4> 

e 


T 

(  Vuj)  A 

-  fie 


-<P 


,1 


-<p 


<P 


2 


Vv^dA 


1  J 


(3.2) 


n  <pds 


-  [  (-V  v20)-u>  ,v2v)  dA 

h  f± 

e 

after  an  integration  by  parts.  Here  n  denotes  the  outward 

pointing  unit  normal  and  (  )  .  indicates  differentiation  with 

o 

respect  to  x..  Since  V  u>  =  0  (see  1.3a),  and  since  from 

2 

(3.1)  it  is  readily  verified  that  V  v  =  0,  it  follows  that  the 
integral  over  on  the  right  hand  of  (3.2)  vanishes. 

Now,  let  us  place  some  concrete  restrictions  on  <p .  Suppose 


that 

(I)  <j>(0)  =  1,  and 


(II) 

tp  =  0  on 

r3’r4’ • *  * 

, T  ( i . e . ,  on  all 

s 

Tj’s  except 

and 

iy. 

Now 

introducing 

the  assumption  a  =  1,2, 

(3.2)  becomes 

l  | 

’  -w,; jV,x  +  U>2V,2" 

A 

4>  = 

•  n  cpds. 

c 

i  JY(e> 

TJl 

-«  -V  2  -  u>  2v  . 

Using  (3.1)  and  the  expansions  (1.4),  (1.5),  we  may  calculate 
directly 


lmt . 
£-*•0 


♦ 

£ 


m 


If  we  assume 
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Q 

(III)  In  the  vicinity  of  0,  V<p  =  0(rp)  where 

g  >  max(0,  — —  — ) ,  then  the  integrand  in  the  first  line  of  the 
a 

definition  (3.2)  of  4*^  is  integrable  over  £),  and  we  have 


(3.3) 


m 


e->-0 


4>  = 

H 

/—> 

3 

> 

"M 

c 

’,1, 

Vv  dA. 


This  is  in  the  form  required  by  (1.1),  with 


X1  =  '<P,1V,1  ‘  <P,2V,2 

X  2  ”  _<P  ,  2V  ,  1  +  ^  ,  lv  ,  2  , 

and  £  =  0,  p  =  0  and  R  =  0. 

In  accord  with  (1.2)  we  may  use  (2.3)  with  w  replacing 

to  to  obtain  an  approximation  km  to  km-  As  usual,  in  practice 

we  would  want  \  to  be  as  smooth  as  possible.  One  strategy  is 

to  let  cp  be  identically  1  in  a  neighborhood  of  0  and  then  let 

II  be  imposed  by  smooth  cut  off  functions.  The  condition  III 

is  automatically  satisfied.  This  strategy  has  the  advantage  of 

eliminating  the  need  to  evaluate  a  singular  integral  near  0 

(since  <p  .  =  <p  _  =  0  there). 

» >  * 

It  may  not  at  first  be  apparent  what  connection  the  express¬ 
ion  (3.3)  has  with  the  energy  release  method.  We  shall  not 
explain  this  in  detail;  however,  let  us  just  note  two  points. 
Firstly,  the  matrix 


i 
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can  be  thought  of  as  a  form  of  stiffness  derivative  with  respect 
to  crack  extension.  (The  function  <p  characterizes  this  crack 
extension.)  While  secondly  the  quantity 


which  appears  in  the  line  integrals  of  (3.2)  can,  by  formally 
putting  v  =  to,  be  modified  to  essentially  give  the  integrand  of 
the  J^-integral  associated  with  (1.3), 


J1(D 


(  1 

r  2  .  2 

“U)  1  +  0)  9 

Jr  2 

-  2gu  _  a)  0 

L  ,1  ,2  J 

•  n  ds . 


I 
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§4.  The  Accuracy  of  the  Approximations  k 

Our  discussion  in  this  section  will  be  brief,  since  much  of 
what  we  said  in  our  earlier  paper  (see  §3.5  of  [  1  ])  applies  in 
the  current  setting  also.  Pursuing  the  same  argument  as  there 
leads,  for  the  techniques  of  both  §2  and  §3,  to  the  error  relation 


(4.1) 


m 


-  k 


m 


(V(io-u)))  V(i|>-iJ>)dA 


where  ^  and  are  the  exact  and  finite  element  solutions  of 

the  auxiliary  problem. 


(4.2) 


-v2t 


V  •  X  +  £  in  n. 


i|i  =  0  on  those  for  which  the  Dirichlet  condition 


^VtJ)  •  n 


(1.3b)  applies. 


A  •  n  +  p  on  those  fj  for  which  the  Neumann 


condition  (1.3c)  applies. 


using  the  notation  of  (1.1)/(1.2).  (In  obtaining  (4.1)  we  have 
implicitly  assumed  that  the  Dirichlet  data  f .  is  exactly 
representable  by  a  finite  element  function.  If  this  is  not  the  case, 
then  the  expression  (4.1)  no  longer  applies.  We  shall  not  discuss 
this  situation  here.)  This  auxiliary  problem  is  of  the  same  form 
as  the  basic  problem  (1.3),  with  only  the  right  hand  sides  of  the 
differential  equation  and  boundary  conditions  changed.  If  we  let 
E(*)  denote  the  strain  energy  expression  associated  with  the  problem 
(1.3)  , 

E( * )  =  f  | V( - ) |2  dA, 

'n 

then  from  (4.1)  we  obtain  the  usual  estimate 

i 
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.  'v  .  ^  1/9  1/9 

(4.3)  k  -k  S  E(u>-w)'l  E(iJ>-^)'l 

1  m  m 

S  min  E( oj-v) ^  min  E(<|>-v*)^^, 
v  v* 

where  v  and  v*  range  over  all  finite  element  functions  which 
vanish  on  those  I\  to  which  the  Dirichlet  condition  (1.3b) 
applies . 

The  conclusion  to  be  drawn  from  (4.3)  is  precisely  the  same 

as  that  which  we  have  already  noted  in  C  1  3 :  the  accuracy  of 

the  post-processed  value  is  related  to  how  well  the  finite  element 

functions  are  able  to  approximate  both  the  solution  w  of  the 

basic  problem  and  the  solution  <J>  of  an  auxiliary  problem. 

The  corner  points  of  ft  usually  give  rise  to  some  form  of 

singular  behavior  in  u  and  \p.  This  seriously  affects  how  well 

these  functions  can  be  approximated  by  the  finite  element  functions. 

Nonetheless,  it  may  be  shown,  in  theory  at  least,  that  for 

suitable  "optimally  refined"  meshes,  with  p*"*1  degree  elements 

min  E(w-v)^^  =  0(N  ) ,  when  N  is  the  number  of  degrees  of 

v 

freedom  of  the  finite  element  model.  The  same  meshes  also  give 

min  E( =  0(N-^^^),  provided  the  right  hand  sides  in  (4.2' 
v* 

are  sufficiently  smooth.  For  such  meshes  then,  we  have  from  (4.3), 
|km~1<m|  =  0(N~^) .  Of  course,  from  a  practical  point  of  view,  the 
constant  multiplying  the  N~^  in  the  0(N~^)  term  is  very 
important.  The  constant  is  related  amongst  other  things,  to  the 
smoothness  of  the  right  hand  sides  of  both  the  basic  problem  (1.3) 
and  the  auxiliary  problem  (4.1).  A  priori,  little  seems  to  be 
able  to  be  said  about  choosing  meshes  that  optimize  its  value. 
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However,  on  a  qualitative  level,  we  again  remark  on  the  importance 
of  choosing  post-processing  procedures  with  "smooth"  X,  £  and 
P’s. 

Let  us  however  remark  on  a  point  that  we  shall  return  to  in 

greater  detail  in  our  next  paper.  Already,  by  concentrating  our 

attention  on  the  estimate  (4.3)  we  have  strayed  somewhat  from  our 

real  objective,  which  is  to  find  high  accuracy  approximations 

to  km-  The  actual  error  in  km  is  given  by  (4.1).  An  optimal 

mesh  strategy,  whether  a  priori  or  adaptive,  should  then  be 

tailored  to  minimizing  the  integral  in  (4.1)  rather  than  the 

right  hand  side  of  the  already  overly  conservative  estimate  (4.3). 

The  estimate  (4.3)  assumes  the  worst  possible  case,  namely  that 

there  is  no  cancellation  in  the  integral  (4.1).  If  there  is 

some  cancellation,  then  k  -  k  will  be  overestimated  by  (4.3). 

mm 
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§  5 .  S^^Nj^erijcal  JExam^les 
§5.1.  The  FEARS  program 

The  calculations  associated  with  the  examples  of  this  section 
were  performed  using  the  FEARS  program.  FEARS  is  a  research 
oriented,  adaptive  finite  element  package  developed  at  the 
University  of  Maryland.  A  detailed  description  of  the  operation 
and  mathematical  background  of  the  program  can  be  found  in  [11]. 

For  the  purposes  of  this  paper,  the  following  few 
remarks  will  suffice.  FEARS  assumes  that  the  region  under  con¬ 
sideration  has  firstly  been  partitioned  into  a  number  of  subregions, 
each  of  which  is  a  curvilinear  quadrilateral.  Within  the  program, 
each  of  these  subregions  is  transformed  by  a  change  of  coordinates 
into  a  unit  square.  The  actual  finite  element  modelling  is  then 
carried  out  on  these  transformed  squares.  Square  bilinear  elements 
are  used.  FEARS  has  an  adaptive  character:  starting  from  an 
initial  coarse  mesh  (usually,  uniform  on  each  of  the  transformed 
squares),  the  program  automatically  selects,  in  a  recursive 
manner,  a  sequence  of  "optimal"  mesh  refinements.  The  program 
allows  the  user  some  freedom  in  choosing  the  criterion  on  which 
the  "optimality"  will  be  based. 

§5.2.  Exam^ljejA :  A^sl  it  ^  j^i^ular  jnembrane 

Let  H  be  the  unit  circle  slit  along  the  positive  axis. 

In  the  notation  of  §1.2  let  be  the  upper  face  of  the  slit, 

r2  the  lower  face  of  the  slit  and  the  circular  portion  of 

the  boundary  of  £).  We  consider  the  following  particular  case 
of  the  model  problem  (1.3): 


(5.1) 


in  0 


on  r. 


on  r, 


on  r. 


(See  Fig.  3).  For  this  problem  the  expansion  (1.4)  about  the 
origin  takes  the  form 

a)  =  k1r1/4sin  |  +  k2r3/4sin  -^  +  k3r5/4sin  0(r7/4) 


(5.2) 


k  -SM/sinC-t  0)\  3k2r-1/Ysln<-7  9>\  Sk.r1'4/31"  * 

Va>  4  \cos(”  0 )/  +  ^  \cos(-i  6)/  4  \cos  i 

+  0(r3/4). 


Using  the  method  of  separation  of  variables,  an  infinite  series 
representation  of  u>  may  be  found.  This  series  can  be  manipulated 
to  give  the  following  exact  values: 


E(u>) 


4.52707, 


kx= -1.35812,  k2  =  .970087,  k3  =  .452707 

(Note  that  the  influence  functions  for  the  k  are  actually 

m 

readily  available  in  the  present  case.  They  are  given  by 


5  =  ir(2m-l)(r 


as  may  be  easily  verified.) 


-(2m-l)~  (2m-l)i 


2ot)sin(  (2m_i)^. 


ki 
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Figure  3. 

Table  1  lists  some  properties  of  a  sequence  of  four  adaptively 
refined  meshes  produced  by  the  FEARS  program  for  the  problem 
(5.1)  with  ft  partitioned  as  in  Figure  3.  FEARS  was  directed  in 
this  instance  to  produce  refinements  that  were  "optimal"  in  the 
strain  energy  sense.  Note  that  the  mesh  I  is  uniform  on  each  of 
the  transformed  squares.  Subsequent  meshes ,  exhibit,  as  would  be 
expected,  quite  severe  refinement  about  the  tip  of  the  slit. 

Using  the  approximate  solutions  corresponding  to  each  of  these 
meshes,  we  calculate  two  sets  of  approximations  lc^,  and  lc3 

to  the  first  three  stress  intensity  factors  of  (5.2).  Firstly, 

< 

we  used  the  Generalized  Influence  Function  Method  of  12  obtained 
by  choosing 
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No.  of  Elements 

m 

MESH  III 

MESH  IV 

(  h  /h  . 

-  max  min 

m 

Subregion  1 

4(1) 

4(1) 

4(1) 

4(1) 

'•  2 

II  tl 

It  II 

It  11 

7(2) 

"  3 

II  II 

II  tl 

It  11 

7(2) 

"  4 

It  II 

It  II 

II  II 

4(1) 

"  5 

II  II 

II  11 

7(2) 

40(2) 

"  6 

II  II 

11  II 

4(1) 

10(2) 

"  7 

II  II 

11  11 

11  II 

13(2) 

"  8 

II  II 

II  II 

11  II 

4(1) 

"  9 

II  II 

10(4) 

13(8) 

100(64) 

"  10 

II  II 

4(1) 

19(4) 

46(64) 

"  11 

II  II 

37(4) 

103(32) 

145(64) 

"  12 

II  II 

25(8) 

40(32) 

133(256) 

Total  No.  of 
Elements 

48 

108 

210 

513 

No.  of  degrees- 
of- freedom 

56 

104 

192 

447 

2 

1  E(  w) 

2.02227 

2.06709 

2.09383 

2.11461 

(31%) 

(24%) 

(18%) 

i _ 

(11%) 

Table  1.  Table  of  mesh  properties  for  Example  A. 


Ci) 


h  ,  h  .  refer  to  the  maximum  and  minimum  mesh  size 

max  min 


on  each  of  the  transformed  squares. 


(ii)  The  quantities  in  parentheses  in  the  last  row  are  the 

_  (  E(u-Z)  \1/2 

"  V  TTu)T“  /  for  the 


relative  energy  norm  errors 
appropriate  meshes. 


(iii)  The  energy  norm  of  the  exact  solution  =  E(w) 


1/2 


2.12769 
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with  v  given  by  (3.1).  The  results  of  these  computations  are 
displayed  in  Table  2. 

The  following  conclusions  may  be  drawn  from  the  results  shown 
in  Tables  1  and  2. 

(a)  Despite  the  presence  of  singularities  at  the  corner 
points  of  ft,  the  sequence  of  adaptively  created  meshes  gives  an 
apparent  rate  of  convergence  for  the  energy  norm  error  in  <*> 
which  is  close  to  the  theoretically  "optimal  refinement"  rate  of 
0(N  where  N  is  the  number  of  degrees-of-freedom  of  the 

finite  element  model.  (The  II  to  III  refinement  step  shows  an 

-.47 

0(N  *  )  rate  and  the  III  to  IV  step  refinement  shows  an 


Mesh 


Influence  Function  Method 


Energy  Release  Method 
(%  relative  error) 


( %  relative  error) 


kl 

k2 

k3 

kl 

k2 

I 

-1.164 

.  9685 

.4529 

-1.170 

.  9698 

.4528 

(31%) 

(14.3%) 

( .16%) 

(.056%) 

(13.8%) 

(  .  025%) 

( .031%) 

II 

-1.243 

.9718 

.4528 

-1.249 

.9734 

(24%) 

(8.5%) 

( . 17%) 

( . 026%) 

(8.0%) 

(  .  34%) 

III 

-1.295 

.9713 

.4532 

-1.300 

.  9726 

.4528 

(18%) 

(4.6%) 

(  .13%) 

(.12%) 

(4.3%) 

( . 26%) 

( . 011%) 

IV 

-1.334 

.9704 

.4533 

-1.338 

.  9700 

.4532 

(11%) 

(1.7%) 

(.036%) 

(  .14%) 

(1.5%) 

(  .012%) 

(  .11%) 

Exact 

-1.358 

.9701 

.4527 

-1.358 

.  9701 

MUlMl 

Value 

i 

■HI 

Table  / 

J.  Table  of  the  computed  values 

of  k  (m=l, 
m 

2,3)  for  Example  A. 

(i)  The  quantity  in  parenthesis  under  each  mesh  label  is 
the  corresponding  relative  energy  norm  error  = 


(ii)  The  exact  values  are  correct  to  the  number  of  significant 
figures  shown. 


_  58 

0(N  *  )  rate).  Had  only  uniform  meshes  been  used  on  the  trans- 

-.125 

formed  squares,  theory  would  have  predicted  an  0(N  *  )  rate. 

In  that  case,  to  obtain  an  accuracy  in  the  energy  norm  comparable 
to  that  of  the  present  mesh  IV,  approxima tely  3  x  10^  degrees- 
of-freedom  would  have  been  required.  This  would  clearly  not  be  a 
practical  proposition . 

(b)  There  seems  to  be  no  significant  difference  between  the 

k  fs  calculated  using  the  Influence  Function  Method  and  those 
m 

calculated  by  the  Energy  Release  Method.  Observe  that  in  both 
cases  appears  to  be  converging  at  a  rate  twice  that  of  the 

energy  norm  error  (i.e.,  at  a  rate  of  approximately  0(N  ^)). 
Notice  also  that  the  coefficients  and  k^  are  very  accurate 

even  for  the  relatively  coarse  meshes  I  and  II.  The  irregular 
be.'.avior  of  the  relative  error  of  k^  and  k^  as  the  meshes 
are  refined  can  probably  be  attributed, at  least  partially,  to 
quadrature  errors  made  in  the  evaluation  of  the  k^.  (We  used 
Gaussian  quadrature).  These  coefficients  are  so  accurate  that 
otherwise  negligible  quadrature  errors  may  now  become  significant 
It  may  be  shown  that  the  auxiliary  functions  \ p  (see  4.2) 
corresponding  to  both  the  Influence  Function  Method  (5.3)  and  the 
Energy  Release  Method  (5.4)  all  take  the  form 

i J)  =  sin(  ( 2m-l  )|~)  (m  =  1,2,3), 

in  the  vicinity  of  the  vertex  0.  Observe  therefore,  that  the 
behavior  of  near  0  becomes  more  "smooth"  as  m  increases. 

This  will  have  consequences  as  regards  the  approximability  of  ij/. 
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Indeed,  the  poorer  accuracy  of  compared  to  k 2  or  kj  can 

be  attributed  to  this  Fact. 


§5.3.  ExajTiple^B :  ^^^i  j^er^r^ane 

2 

Let  ft  be  the  square  (-1,1)  which  has  been  slit  along  the 
axis.  In  the  notation  of  §1.2  let  be  the  upper  face  of 

the  slit,  r2  the  lower  face  of  the  slit  and  let  r3,...,r2  be 
some  labelling  of  the  remaining  straight  line  segments  making  up 
the  boundary  of  ft.  We  consider  the  following  particular  case  of 
(1.3)  : 


(5.3) 


0  on  ft 

o  on  r1 

o  on  r2 

x2  and  r  ^ , . . . j  F  y ■ 


(See  Fig.  4).  For  this  problem  the  expansions  (5.2)  remain  valid 
about  the  origin.  Using  the  technique  of  conformal  mappings  and 
the  method  of  separation  of  variables,  an  infinite  series  repre¬ 
sentation  of  w  may  be  found.  From  this  series  it  can  be  deduced 
that 

k  =  -.353768,  k?  =  .735947,  k3  =  .5608491. 

Unfortunately,  this  series  does  not  provide  an  effective  means  for 
calculating  the  exact  energy  E(u»)  of  the  exact  solution  u>. 
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Figure  4. 

Table  3  lists  some  of  the  characteristics  of  a  sequence  of 
four  adaptively  refined  meshes  which  were  created  by  FEARS  for 
this  problem.  The  initial  partitioning  of  0,  into  subregions 
is  illustrated  in  Figure  4.  FEARS  was  directed  in  this  instance 
to  produce  refinements  that  were  "optimal"  in  the  strain  energy 
sense.  Again  the  first  mesh  is  uniform  on  each  of  the  transformed 
squares.  The  subsequent  meshes  show,  as  expected,  refinement 
around  the  slit  tip.  Using  the  approximate  solutions  correspond¬ 
ing  to  each  of  these  meshes,  we  again  calculate  two  sets  of 
approximations  k^,  Jc2 ,  ]<3  to  the  first  three  stress  intensity 
factors  of  (5.2).  Firstly,  we  employ  the  Generalized  Influence 
Method  of  §2.  In  particular  we  choose 


No.  of  elements 

1  MESH  I 

MESH  II 

MESH  III 

MESH  IV 

Ch  /h  .  ) 
max   mm 

Subregion  1 

4(1) 

4(1) 

4(1) 

4(1) 

It 

2 

II  II 

II  tl 

II  11 

11  11 

II 

3 

II  II 

II  It 

II  It 

II  II 

11 

4 

II  II 

II  II 

II  11 

II  11 

tl 

5 

If  If 

II  II 

II  II 

11  It 

II 

6 

It  II 

11  11 

67(8) 

130(64) 

II 

7 

II  II 

11  11 

19(4) 

19(4) 

II 

8 

II  II 

II  II 

4(1) 

4(1) 

f? 

9 

II  II 

11  II 

4(1) 

4(1) 

II 

10 

II  II 

67(8) 

70(8) 

91(64) 

II 

11 

II  II 

34(4) 

43(16) 

76(16) 

II 

12 

II  II 

37(4) 

40(4) 

46(4) 

II 

13 

II  II 

4(1) 

4(1) 

4(1) 

II 

14 

It  II 

II  11 

II  11 

11  11 

II 

15 

l 

II  II 

II  11 

II  11 

7(2) 

II 

16 

II  11 

tl  It 

1!  11 

7(2) 

Total  No 
Elements 

.  of 

1 

64 

190 

283 

412 

No.  of  degrees- 
of-freedom 

48 

159 

240 

348 

E(u) 

1/2  j 

1.93212 

1.92293 

1.92038  | 

1.91888 

Table  3.  Table  of  mesh  properties  for  Example  B. 

,  h  .  ref ex  to  the  maximum  and  minimum  mesh 

max  min 


size  on  each  of  the  transformed  squares. 
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where 


X(x1)X(x2) 


ttC  2m- 1) 


sin 


X(t) 


|t| 

|t| 


< 


> 


1 

2 

1 

2 


i 

i 

i 


in  (2.2).  Notice  that  in  contrast  to  the  corresponding  choice  in 
Example  A,  we  here  had  to  require  that  tp  =  0  on  the  non-slit 
part  of  3S).  We  have  done  this  by  using  the  cut  off  function 
X(x^)X(x2>.  The  expression  used  to  evaluate  is  the  particular 

form  of  (2.2) 


k 

m 


•  ft 


V  cpiodA 


l  Vip  •  nx.ds  . 

j  =  3  Jr.  l 

J  3 


The  second  set  of  approximations  is  based  on  the  Generalized 
Energy  Release  Method  of  §3  with  the  function  <p  of  (3.3) 
given  by 


where 


<p(x1,x2)  =  cp(x1)(p(x2)  , 


(  1  0  <  1 1 1  <  | 
\  1  “  2<|t|  -|)  \  5  1 1 1  <  1 


The  results  of  these  computations  are  displayed  in  Table  4. 

Just  as  in  Example  A,  we  see  that  for  this  example  there  are 
no  significant  differences  in  accuracy  between  the  influence  func¬ 
tion  and  energy  release  methods.  We  also  see  here  the  high  accuracy 

of  k2  and  k^  even  for  the  coarse  meshes  I  and  II.  Note  also  that 

~  -1 

the  convergence  of  k^  is  consistent  with  the  0(N  )  rate  that  we 

would  have  for  a  theoretically  optimal  refinement. 


29 


Mesh 

Influence  Function 

(%  relative  err 

Method 

*or) 

Energy  Release  Method 

(%  relative  error) 

I 

V 

K1 

*2 

£ 

k3 

kl 

k2 

k3 

-.4309 

(22%) 

.7400 

(  .  56%) 

.  5602 

( .12%) 

-.4203 

(19%) 

.  7361 

(  .  024%) 

.  5627 

( .33%) 

II 

-.  3941 

(11%) 

.  7364 

(  .  067%) 

.  5605 

(  .  053%) 

-.3924 

(11%) 

.  7357 

( . 032%) 

.  5610 

( .029%) 

III 

-.3750 

(6.0%) 

.7364 

(  .068%) 

.  5607 

( .033%) 

-.3748 

(6.0%) 

.  7361 

( .020%) 

IV 

1 

-.3671 

(3.8%) 

.7363 

(  .052%) 

.  5607 

( .027%) 

-.3669 

(3.7%) 

.  7359 

(  .0071%) 

.5610 

( .027%) 

Exact 

Value 

-.3538 

.  7359 

.  5608 

-.  3538 

.  7359 

.  5608 

Table  4.  Table  of  computed  values  of  km  Cm  =  1,2,3)  for  Example  B. 


The  exact  values  are  correct  to  the  number  of  significant 
figures  shown. 
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§5.4.  Conclusions 

For  these  examples  we  see  that  the  use  of  post-processing 
techniques  allied  with  an  appropriate  adaptive  mesh  selection 
algorithm,  makes  it  possible  to  obtain  the  leading  stress  intensity 
factor  with  say  a  five-percent  accuracy,  while  using  only  a  moderate 
number  of  degrees-of-freedom  in  the  finite  element  model  (around 
250  for  our  examples).  This ,  even  though  the  strain  energy  norm 
error  for  such  a  mesh  may  still  be  of  the  order  of  20%  (as  in 
Example  A) .  This  latter  observation  also  shows  the  importance  of 
being  specific  when  speaking  of  the  accuracy  of  a  finite  element 
solution  -  a  solution  giving  acceptable  accuracy  for  one  quantity, 
may  not  be  acceptable  for  another. 
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